The present invention is directed to the field of signal processing, and, more particularly, is directed to systems and methods for signal averaging.
Digital signal processing techniques are frequently employed to enhance a desired signal in a wide variety of applications, such as health care, communications and avionics, to name a few. Signal enhancement includes smoothing, filtering and prediction. These processing techniques each operate on a block of input signal values in order to estimate the signal at a specific point in time. FIG. 1 illustrates that smoothing, filtering and prediction can be distinguished by the time at which an output value is generated relative to input values. Shown in FIG. 1 is a time axis 100 and a block 101 of input signal values depicted in this example as occurring within a time window between points tmin and tmax. Specifically, the block 101 includes a set of discrete input values {vi; i=1, 2, . . . n} occurring at a corresponding set of time points {ti; i=1, 2, . . . n}. A smoother operates on the block 101 of input values to estimate the signal at a time point, ts 102 between tmin and tmax. That is, a smoother generates an output value based upon input values occurring before and after the output value. A filter operates on the block 101 of input values to estimate the signal at a time tf 104, corresponding to the most recently occurring input value in the block 101. That is, a filter generates a forward filtered output value at the time tf based upon input values occurring at, and immediately before, the output value. A filter also operates on the block 101 to estimate the signal at a time tb 105 at the beginning of the block 101 to generate a backward filtered value. A forward predictor operates on the block of input values 101 to estimate the signal at time tpf 106, which is beyond the most recently occurring value in the block 101. That is, a forward predictor generates a forward predicted output value based upon input values occurring prior to the output value. A backward predictor operates on the block 101 of input values to estimate the signal at time tpb 108, which is before the earliest occurring value in the block 101. That is, a backward predictor generates a backward predicted output value based upon input values occurring after the output value.
A common smoothing technique uses an average to fit a constant, vA, to a set of data values, {vi; i=1, 2, . . . , n}:                               v          A                =                              1            n                    ·                                    ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                          v              i                                                          (        1        )            
A generalized form of equation (1) is the weighted average                               v          WA                =                                            ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                          w                i                            ·                              v                i                                                                        ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                          w              i                                                          (        2        )            
Here, each value, vi, is scaled by a weight, wi, before averaging. This allows data values to be emphasized and de-emphasized relative to each other. If the data relates to an input signal, for example, values occurring during periods of low signal confidence can be given a lower weight and values occurring during periods of high signal confidence can be given a higher weight.
FIG. 2A illustrates the output of a constant mode averager, which utilizes the weighted average of equation (2) to process a discrete input signal, {vi; i an integer} 110. The input signal 110 may be, for example, a desired signal corrupted by noise or a signal having superfluous features. The constant mode averager suppresses the noise and unwanted features, as described with respect to FIG. 5, below. A first time-window 132 defines a first set, {vi; i=1, 2, . . . , n}, of signal values, which are averaged together to produce a first output value, z1 122. A second time-window 134, shifted from the previous window 132, defines a second set {vi; i=2, 3, . . . , n+1}of signal values, which are also averaged together to produce a second output value z2 124. In this manner, a discrete output signal, {zj; j an integer} 120 is generated from a moving weighted average of a discrete input signal {vi; i an integer} 110, where:                               z          j                =                              ∑                          i              =              j                                      n              +              j              -              1                                ⁢                      xe2x80x83                    ⁢                                    w              i                        ⁢                                          v                i                            /                                                ∑                                      i                    =                    j                                                        n                    +                    j                    -                    1                                                  ⁢                                  xe2x80x83                                ⁢                                  w                  i                                                                                        (        3        )            
A common filtering technique computes a linear fit to a set of data values, {vi; i=1, 2, . . . , n}:
{circumflex over (v)}i=xcex1xc2x7ti+xcex2xe2x80x83xe2x80x83(4)
where xcex1 and xcex2 are constants and ti is the time of occurrence of the ith value. FIG. 2B illustrates the output of a linear mode averager, which uses the linear fit of equation (4) to process a discrete input signal, {vi; i an integer} 110. The input signal 110 may be, for example, a desired signal with important features corrupted by noise. The linear mode averager reduces the noise but tracks the important features, as described with respect to FIG. 6 below. A first time-window 132 defines a first set, {vi; i=1, 2, . . . , n}, of signal values. A linear fit to these n values is a first line 240, and the value along this line at max {t1, t2, . . . , tn} is equal to a first output value, z1 222. A second time-window 134 shifted from the previous window 132 defines a second set, {vi; i=2, 3, . . . , n+1 }, of signal values. A linear fit to these n values is a second line 250, and the value along this line at max {t2, t3, . . . , tn+1} is equal to a second output value, z2 224. In this manner, a discrete output signal, {zj; j an integer} 220 is generated from a moving linear fit of a discrete input signal {vi; i an integer}, where:                               z          j                =                                            α              j                        ·                          t                              n                +                j                -                1                            MAX                                +                      β            j                                              (5a)                                          t                      n            +            j            -            1                    MAX                =                  max          ⁢                      {                                          t                j                            ,                              t                                  j                  +                  1                                            ,              …              ⁢                              xe2x80x83                            ,                              t                                  n                  +                  j                  -                  1                                                      }                                              (5b)            
In general, the time windows shown in FIGS. 2A-2B may be shifted from each other by more than one input value, and values within each time window may be skipped, i.e., not included in the average. Further, the ti""s may not be in increasing or decreasing order or uniformly distributed, and successive time windows may be of different sizes. Also, although the discussion herein refers to signal values as the dependent variable and to time as the independent variable to facilitate disclosure of the present invention, the concepts involved are equally applicable where the variables are other than signal values and time. For example, an independent variable could be a spatial dimension and a dependent variable could be an image value.
The linear mode averager described with respect to FIG. 2B can utilize a xe2x80x9cbestxe2x80x9d linear fit to the input signal, calculated by minimizing the mean-squared error between the linear fit and the input signal. A weighted mean-squared error can be described utilizing equation (4) as:                               ϵ          ⁡                      (                          α              ,              β                        )                          =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                                                      w                  i                                ⁡                                  (                                                            v                      i                                        -                                                                  v                        ^                                            i                                                        )                                            2                        /                                          ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                              w                i                                                                        (6a)                                          ϵ          ⁡                      (                          α              ,              β                        )                          =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                                                      w                  i                                ⁡                                  [                                                            v                      i                                        -                                          (                                                                        α                          ·                                                      t                            i                                                                          +                        β                                            )                                                        ]                                            2                        /                                          ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                              w                i                                                                        (6b)            
Conventionally, the least-mean-squared (LMS) error is calculated by setting the partial derivatives of equation (6b) with respect to xcex1 and xcex2 to zero:                                           ∂                          ∂              α                                ⁢                      xe2x80x83                    ⁢                      ϵ            ⁡                          (                              α                ,                β                            )                                      =        0                            (7a)                                                      ∂                          ∂              β                                ⁢                      xe2x80x83                    ⁢                      ϵ            ⁡                          (                              α                ,                β                            )                                      =        0                            (7b)            
Substituting equation (6b) into equation (7b) and taking the derivative yields:                                           -            2                    ⁢                                    ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                                            w                  i                                ⁡                                  [                                                            v                      i                                        -                                          (                                                                        α                          ·                                                      t                            i                                                                          +                        β                                            )                                                        ]                                            /                                                ∑                                      i                    =                    1                                    n                                ⁢                                  xe2x80x83                                ⁢                                  w                  i                                                                    =        0                            (8)            
Solving equation (8) for xcex2 and substituting the expression of equation (2) yields:                     β        =                                                            ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                                                w                  i                                ·                                  v                  i                                                                                    ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                              w                i                                              -                      α            ⁢                                                            ∑                                      i                    =                    1                                    n                                ⁢                                  xe2x80x83                                ⁢                                                      w                    i                                    ·                                      t                    i                                                                                                ∑                                      i                    =                    1                                    n                                ⁢                                  xe2x80x83                                ⁢                                  w                  i                                                                                        (9a)            
xcex2=xcexdWAxe2x88x92xcex1xc2x7tWAxe2x80x83xe2x80x83(9b)
where the weighted average time, tWA, is defined as:                               t          WA                =                                            ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                          w                i                            ·                              t                i                                                                        ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                          w              i                                                          (        10        )            
Substituting equation (9b) into equation (4) gives:
{circumflex over (v)}i=xcex1(tixe2x88x92tWA)+vWAxe2x80x83xe2x80x83(11)
Substituting equation (11) into equation (6a) and rearranging terms results in:                               ϵ          ⁡                      (                          α              ,              β                        )                          =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                                                      w                  i                                ⁡                                  [                                                            (                                                                        v                          i                                                -                                                  v                          WA                                                                    )                                        -                                          α                      ·                                              (                                                                              t                            i                                                    -                                                      t                            WA                                                                          )                                                                              ]                                            2                        /                                          ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                              w                i                                                                        (        12        )            
Changing variables in equation (12) gives:                               ϵ          ⁡                      (                          α              ,              β                        )                          =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                                                      w                  i                                ⁡                                  (                                                            v                      i                      xe2x80x2                                        -                                          α                      ·                                              t                        i                        xe2x80x2                                                                              )                                            2                        /                                          ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                              w                i                                                                        (        13        )            
where:
vxe2x80x2i=vixe2x88x92vWAxe2x80x83xe2x80x83(14a)
txe2x80x2i=tixe2x88x92tWAxe2x80x83xe2x80x83(14b)
Substituting equation (13) into equation (7a) and taking the derivative yields                                           -            2                    ⁢                                    ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                          w                i                            ⁢                                                                    t                    i                    xe2x80x2                                    ⁡                                      (                                                                  v                        i                        xe2x80x2                                            -                                              α                        ·                                                  t                          i                          xe2x80x2                                                                                      )                                                  /                                                      ∑                                          i                      =                      1                                        n                                    ⁢                                      xe2x80x83                                    ⁢                                      w                    i                                                                                      =        0                            (        15        )            
Solving equation (15) for xcex1 gives:                     α        =                                            ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                          w                i                            ⁢                              v                i                xe2x80x2                            ⁢                                                t                  i                  xe2x80x2                                /                                                      ∑                                          i                      =                      1                                        n                                    ⁢                                      xe2x80x83                                    ⁢                                      w                    i                                                                                                          ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                          w                i                            ⁢                                                t                  i                                      xe2x80x2                    2                                                  /                                                      ∑                                          i                      =                      1                                        n                                    ⁢                                      xe2x80x83                                    ⁢                                      w                    i                                                                                                          (        16        )            
Substituting equations (14a, b) into equation (16) results in:                     α        =                                            ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                                            w                  i                                ⁡                                  (                                                            v                      i                                        -                                          v                      WA                                                        )                                            ⁢                                                (                                                            t                      i                                        -                                          t                      WA                                                        )                                /                                                      ∑                                          i                      =                      1                                        n                                    ⁢                                      xe2x80x83                                    ⁢                                      w                    i                                                                                                          ∑                              i                =                1                            n                        ⁢                          xe2x80x83                        ⁢                                                                                w                    i                                    ⁡                                      (                                                                  t                        i                                            -                                              t                        WA                                                              )                                                  2                            /                                                ∑                                      i                    =                    1                                    n                                ⁢                                  xe2x80x83                                ⁢                                  w                  i                                                                                        (17a)                                α        =                              σ            vt            2                                σ            tt            2                                              (17b)            
where:                               σ          vt          2                =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                                    w                i                            ⁡                              (                                                      v                    i                                    -                                      v                    WA                                                  )                                      ⁢                                          (                                                      t                    i                                    -                                      t                    WA                                                  )                            /                                                ∑                                      i                    =                    1                                    n                                ⁢                                  xe2x80x83                                ⁢                                  w                  i                                                                                        (18a)                                          σ          tt          2                =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                                                      w                  i                                ⁡                                  (                                                            t                      i                                        -                                          t                      WA                                                        )                                            2                        /                                          ∑                                  i                  =                  1                                n                            ⁢                              xe2x80x83                            ⁢                              w                i                                                                        (18b)            
Finally, substituting equation (17b) into equation (11) provides the equation for the least-mean-square (LMS) linear fit to {vi; i=1, 2, . . . , n}:                                           v            ^                    i                =                                                            σ                vt                2                                            σ                tt                2                                      ⁢                          xe2x80x83                        ⁢                          (                                                t                  i                                -                                  t                  WA                                            )                                +                      v            WA                                              (        19        )            
FIG. 3 provides one comparison between the constant mode averager, described above with respect to FIG. 2A and equation (2), and the linear mode averager, described above with respect to FIG. 2B and equation (19). Shown in FIG. 3 are input signal values {vi; i=1, 2, . . . , n } 310. The constant mode averager calculates a constant 320 for these values 310, which is equal to vWA, the weighted average of the input values vi. Thus, the constant mode averager output 340 has a value vWA. For comparison to the linear mode averager, the constant mode averager output can be conceptualized as an estimate of the input values 310 along a linear fit 350, evaluated at time tWA. The linear mode averager may be thought of as calculating a LMS linear fit, {circumflex over (v)}i 330 to the input signal values, vi 310. The linear mode averager output 350 has a value, vWLA. The linear mode averager output is an estimate of the input values 310 along the linear fit 330, described by equation (19), evaluated at an index i such that ti=tMAX:                               v          WLA                =                                                            σ                vt                2                                            σ                tt                2                                      ⁢                          xe2x80x83                        ⁢                          (                                                t                  MAX                                -                                  t                  WA                                            )                                +                      v            WA                                              (        20        )            
where:
tMAX=max{t1, t2, . . . , tn}xe2x80x83xe2x80x83(21)
As illustrated by FIG. 3, unlike the constant mode averager, the linear mode averager is sensitive to the input signal trend. That is, the constant mode averager provides a constant fit to the input values, whereas the linear mode averager provides a linear fit to the input values that corresponds to the input value trend. As a result, the output of the linear mode averager output responds faster to changes in the input signal than does the output of the constant mode averager. The time lag or delay between the output of the constant mode averager and the output of the linear mode averager can be visualized by comparing the time difference 360 between the constant mode averager output value 340 and the linear mode averager output value 350.
FIGS. 4-6 illustrate further comparisons between the constant mode averager and the linear mode averager. FIG. 4 depicts a noise-corrupted input signal 410, which increases in frequency with time. FIGS. 5-6 depict the corresponding noise-free signal 400. FIG. 5 also depicts the constant mode averager output 500 in response to the input signal 410, with the noise-free signal 400 shown for reference. FIG. 6 depicts the linear mode averager output 600 in response to the input signal 410, with the noise-free signal 400 also shown for reference. As shown in FIG. 5, the constant mode averager output 500 suppresses noise from the input signal 410 (FIG. 4) but displays increasing time lag and amplitude deviation from the input signal 400 as frequency increases. As shown in FIG. 6, the linear mode averager output 600 tends to track the input signal 400 but also tracks a portion of the noise on the input signal 410.
FIGS. 4-6 suggest that it would be advantageous to have an averager that has variable characteristics between those of the linear mode averager and those of the constant mode averager, depending on signal confidence. Specifically, it would be advantageous to have a variable mode averager that can be adjusted to track input signal features with a minimal output time lag when signal confidence is high and yet adjusted to smooth an input signal when signal confidence is low. Further, it would be advantageous to have a variable mode averager that can be adjusted so as not to track superfluous input signal features regardless of signal confidence.
One aspect of the present invention is a variable mode averager having a buffer that stores weighted input values. A mode input specifies a time value relative to the input values. A processor is coupled to the buffer, and the processor is configured to provide an estimate of the input values that corresponds to the time value. In a particular embodiment, the mode input is adjustable so that the estimate varies between that of a smoother and that of a forward predictor of the input values. In another embodiment, the mode input is adjustable so that the estimate varies between that of a smoother and that of a filter of the input values. In yet another embodiment, the mode input is adjustable so that the estimate varies between that of an average of the input values and that of a filter of the input values. The mode input may be adjustable based upon a characteristic associated with the input values, such as a confidence level. In one variation of that embodiment, the estimate can be that of a smoother when the confidence level is low and that of a filter when the confidence level is high. The estimate may occur along a curve-fit of the input values at the time value. In one embodiment, the curve-fit is a linear LMS fit to the input values.
Another aspect of the present invention is a signal averaging method. The method includes identifying signal values and determining weights corresponding to the signal values. The method also includes computing a trend of the signal values adjusted by the weights. Further, the method includes specifying a time value relative to the signal values based upon a characteristic associated with the signal values and estimating the signal values based upon the trend evaluated at the time value. The method may also incorporate the steps of determining a confidence level associated with the signal values and specifying the time value based upon the confidence level. In one embodiment, the trend is a linear LMS fit to the signal values adjusted by the weights. In that case, the time value may generally correspond to the maximum time of the signal values when the confidence level is high and generally correspond to the weighted average time of the signal values when the confidence level is low.
Yet another aspect of the present invention is a signal averaging method having the steps of providing an input signal, setting a mode between a first mode value and a second mode value and generating an output signal from an estimate of the input signal as a function of said mode. The output signal generally smoothes the input signal when the mode is proximate the first mode value, and the output signal generally tracks the input signal when the mode is proximate the second mode value. The method may also include determining a characteristic of the input signal, where the setting step is a function of the characteristic. In one embodiment, the characteristic is a confidence level relating to the input signal. In another embodiment, the setting step incorporates the substeps of setting the mode proximate the first mode value when the confidence level is low and setting the mode proximate the second mode value when the confidence level is high. In another embodiment, the input signal is a physiological measurement and the setting step comprises setting the mode proximate the first mode value when the measurement is corrupted with noise or signal artifacts and otherwise setting the mode proximate the second mode value so that the output signal has a fast response to physiological events.
A further aspect of the present invention is a signal averager having an input means for storing signal values, an adjustment means for modifying the signal values with corresponding weights, a curve fitting means for determining a trend of the signal values, and an estimate means for generating an output value along the trend. The signal averager may further have a mode means coupled to the estimate means for variably determining a time value at which to generate the output value.